1. 第二章知识补充: 多项式回归
我们在前面讲的都是一般线性回归,即使用的假设函数是一元一次方程,也就是二维平面上的一条直线。
但是很多时候可能会遇到直线方程无法很好的拟合数据的情况,这个时候可以尝试使用多项式回归。
多项式回归中,加入了特征的更高次方(例如平方项或立方项),也相当于增加了模型的自由度,用来捕获数据中非线性的变化。添加高阶项的时候,也增加了模型的复杂度。随着模型复杂度的升高,模型的容量以及拟合数据的能力增加,可以进一步降低训练误差,但导致过拟合的风险也随之增加(后面会专门讨论出现过拟合的情况)。
2. 1 多项式回归的一般形式
在多项式回归中,最重要的参数是最高次方的次数。设最高次方的次数为n,且只有一个特征时,其多项式回归的方程为:
如果令x_0=1x0=1,在多样本的情况下,可以写成向量化的形式:
其中𝑋是大小为m⋅(n+1)的矩阵,θ是大小为(n+1)⋅1的矩阵。
在这里虽然只有一个特征x以及x的不同次方,但是也可以将x的高次方当做一个新特征。与多元回归分析唯一不同的是,这些特征之间是高度相关的,而不是通常要求的那样是相互对立的。
在这里有个问题在刚开始学习线性回归的时候困扰了自己很久:如果假设中出现了高阶项,那么这个模型还是线性模型吗?此时看待问题的角度不同,得到的结果也不同。如果把上面的假设看成是特征xx的方程,那么该方程就是非线性方程;如果看成是参数𝜃θ的方程,那么xx的高阶项都可以看做是对应𝜃θ的参数,那么该方程就是线性方程。很明显,在线性回归中采用了后一种解释方式。因此多项式回归仍然是参数的线性模型。
3. 2 多项式回归的实现
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
下是使用的数据是使用 y=x2+2 并加入一些随机误差生成的,只取了10个数据点:
# 构造数据,数据可视化展示
data = np.array([[ -2.95507616, 10.94533252],
[ -0.44226119, 2.96705822],
[ -2.13294087, 6.57336839],
[ 1.84990823, 5.44244467],
[ 0.35139795, 2.83533936],
[ -1.77443098, 5.6800407 ],
[ -1.8657203 , 6.34470814],
[ 1.61526823, 4.77833358],
[ -2.38043687, 8.51887713],
[ -1.40513866, 4.18262786]])
X = data[:, 0].reshape(-1, 1) # 将array转换成矩阵
y = data[:, 1].reshape(-1, 1)
plt.plot(X, y, "b.")
plt.xlabel('X')
plt.ylabel('y')
plt.show()
3.1. 2.1 直线方程的拟合
下面先用直线方程拟合上面的数据点
lin_reg = LinearRegression()
lin_reg.fit(X, y)
print(lin_reg.intercept_, lin_reg.coef_) # [ 4.97857827] [[-0.92810463]]
X_plot = np.linspace(-3, 3, 1000).reshape(-1, 1)
# 可以使用两种方法用于模型预测
# y_plot = np.dot(X_plot, lin_reg.coef_.T) + lin_reg.intercept_
y_plot = lin_reg.predict(X_plot)
plt.plot(X_plot, y_plot,"red")
plt.plot(X, y, 'b.')
plt.xlabel('X')
plt.ylabel('y')
# 使用mse衡量其误差值:
y_pre = lin_reg.predict(X)
mean_squared_error(y, y_pre)
# 3.3363076332788495
3.2. 2.2 使用多项式方程
为了拟合2次方程,需要有特征x2的数据,这里可以使用函数"PolynomialFeatures"来获得:
sklearn 的 PolynomialFeatures 的用法
使用 sklearn.preprocessing.PolynomialFeatures
这个类可以进行特征的构造,构造的方式就是特征与特征相乘(自己与自己,自己与其他人),这种方式叫做使用多项式的方式。 例如:有 𝑎、𝑏 两个特征,那么它的 2 次多项式的次数为 1,a,b,a2,ab,b2。
PolynomialFeatures
这个类有 3 个参数:
- degree:控制多项式的次数;
- interaction_only:默认为 False,如果指定为 True,那么就不会有特征自己和自己结合的项,组合的特征中没有 a2 和 b2 ;
- include_bias:默认为 True 。如果为 True 的话,那么结果中就会有 0 次幂项,即全为 1 这一列。
X = np.arange(6).reshape(3, 2)
X
# 输出结果
array([[0, 1],
[2, 3],
[4, 5]])
from sklearn.preprocessing import PolynomialFeatures
# 设置多项式阶数为2,其他值默认
# degree 多项式阶数
poly = PolynomialFeatures(degree=2)
res = poly.fit_transform(X)
res
# 输出结果
array([[ 1., 0., 1., 0., 0., 1.],
[ 1., 2., 3., 4., 6., 9.],
[ 1., 4., 5., 16., 20., 25.]])
使用函数"PolynomialFeatures"获取二次方项:
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
print(X_poly)
# 输出结果
[[-2.95507616 8.73247511]
[-0.44226119 0.19559496]
[-2.13294087 4.54943675]
[ 1.84990823 3.42216046]
[ 0.35139795 0.12348052]
[-1.77443098 3.1486053 ]
[-1.8657203 3.48091224]
[ 1.61526823 2.60909145]
[-2.38043687 5.66647969]
[-1.40513866 1.97441465]]
利用上面的数据做线性回归分析:
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
print(lin_reg.intercept_, lin_reg.coef_)
# [ 2.60996757] [[-0.12759678 0.9144504 ]]
X_plot = np.linspace(-3, 3, 1000).reshape(-1, 1)
X_plot_poly = poly_features.fit_transform(X_plot)
y_plot = lin_reg.predict(X_plot_poly)
plt.plot(X_plot, y_plot, 'red')
plt.plot(X, y, 'b.')
plt.show()
第3行得到了训练后的参数,即多项式方程为h=-0.13x+0.91x^2+2.61h=−0.13x+0.91x2+2.61(结果中系数的顺序与𝑋中特征的顺序一致),如下图所示:
# 使用mse衡量其误差值:
y_pre = lin_reg.predict(X_poly)
mean_squared_error(y, y_pre)
# 0.07128562789085331
利用多项式回归,代价函数MSE的值下降到了0.07。
通过观察代码,可以发现训练多项式方程与直线方程唯一的差别是输入的训练集𝑋的差别。在训练直线方程时直接输入了𝑋的值,在训练多项式方程的时候,还添加了我们计算出来的x2这个“新特征”的值(由于x2完全是由x的值确定的,因此严格意义上来讲此时该模型只有一个特征x)。
此时有个非常有趣的问题:假如一开始得到的数据就是上面代码中"X_poly"的样子,且不知道x1与x2之间的关系。此时相当于我们有10个样本,每个样本具有x1,x2两个不同的特征。这时假设函数为:
直接按照二元线性回归方程来训练,也可以得到上面同样的结果(𝜃的值)。如果在相同情况下,收集到了新的数据,可以直接带入上面的方程进行预测。
唯一不同的是,我们不知道x2=x^2_1x2=x12这个隐含在数据内部的关系,所有也就无法画出上图中的这条曲线。一旦了解到了这两个特征之间的关系,数据的维度就从3维下降到了2维(包含截距项\theta_0θ0)。
4. 3 持续降低训练误差与过拟合
在上面实现多项式回归的过程中,通过引入高阶项x2,训练误差从3.34下降到了0.07,减小了将近50倍。那么训练误差是否还有进一步下降的空间呢?
- 答案是肯定的,通过继续增加更高阶的项,训练误差可以进一步降低。通过尝试,当最高阶项为x12时,训练误差为3.87e-23,几乎等于0了。
下面是测试不同degree的过程:
# 定义模型训练函数
def try_degree(degree, X, y):
poly_features_d = PolynomialFeatures(degree=degree, include_bias=False)
X_poly_d = poly_features_d.fit_transform(X)
lin_reg_d = LinearRegression()
lin_reg_d.fit(X_poly_d, y)
return {'X_poly': X_poly_d, 'intercept': lin_reg_d.intercept_, 'coef': lin_reg_d.coef_}
degree2loss_paras = []
for i in range(2, 20):
paras = try_degree(i, X, y)
# 自己实现预测值的求解
h = np.dot(paras['X_poly'], paras['coef'].T) + paras['intercept']
_loss = mean_squared_error(h, y)
degree2loss_paras.append({'d': i, 'loss': _loss, 'coef': paras['coef'], 'intercept': paras['intercept']})
查看最小模型参数:
min_index = np.argmin(np.array([i['loss'] for i in degree2loss_paras]))
min_loss_para = degree2loss_paras[min_index]
print(min_loss_para)
# 输出结果
{'d': 12,
'loss': 3.8764202841976227e-23,
'coef': array([[ 1.17159189, 8.60674192, -4.91798703, -4.18378115, 3.79426131, -8.56026107, -6.94465715, 5.03891035, 4.08870088, -0.30369348, -0.6635493 , -0.11314395]]),
'intercept': array([1.63695924])}
对最小模型可视化展示:
X_plot = np.linspace(-3, 1.9, 1000).reshape(-1, 1)
poly_features_d = PolynomialFeatures(degree=min_loss_para['d'], include_bias=False)
X_plot_poly = poly_features_d.fit_transform(X_plot)
y_plot = np.dot(X_plot_poly, min_loss_para['coef'].T) + min_loss_para['intercept']
plt.plot(X_plot, y_plot, 'red', label="degree12")
plt.plot(X, y, 'b.', label="X")
plt.legend(loc='best')
plt.show()
此时函数图像穿过了每一个样本点,所有的训练样本都落在了拟合的曲线上,训练误差接近与0。 可以说是近乎完美的模型了。但是,这样的曲线与我们最开始数据的来源(一个二次方程加上一些随机误差)差异非常大。
如果从相同来源再取一些样本点,使用该模型预测会出现非常大的误差。类似这种训练误差非常小,但是新数据点的测试误差非常大的情况,就叫做模型的过拟合。过拟合出现时,表示模型过于复杂,过多考虑了当前样本的特殊情况以及噪音(模型学习到了当前训练样本非全局的特性),使得模型的泛化能力下降。
防止模型过拟合是机器学习领域里最重要的问题之一。鉴于该问题的普遍性和重要性,在满足要求的情况下,能选择简单模型时应该尽量选择简单的模型。